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Abstract. A fourth-order time-splitting Laguerre-Hermite pseudospectral method is introduced 
for Bose-Einstein condensates (BEC) in 3-D with cylindrical symmetry. The method is explicit, un- 
conditionally stable, time reversible and time transverse invariant. It conserves the position density, 
and is spectral accurate in space and fourth-order accurate in time. Moreover, the new method 
has two other important advantages: (i) it reduces a 3-D problem with cylindrical symmetry to an 
effective 2-D problem; (ii) it solves the problem in the whole space instead of in a truncated artificial 
computational domain. The method is applied to vector Gross-Pitaevskii equations (VGPEs) for 
multi-component BECs. Extensive numerical tests are presented for 1-D GPE, 2-D GPE with radial 
symmetry, 3-D GPE with cylindrical symmetry as well as 3-D VGPEs for two-component BECs to 
show the efficiency and accuracy of the new numerical method. 
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1. Introduction. Since its realization in dilute bosonic atomic gases E], 
Bose-Einstein condensation of alkali atoms and hydrogen has been produced and 
studied extensively in the laboratory |26j . and has spurred great excitement in the 
atomic physics community and renewed the interest in studying the collective dy- 
namics of macroscopic ensembles of atoms occupying the same one-particle quantum 
state [23 El I2H • Theoretical predictions of the properties of a BEC like the density 
profile collective excitations |2J| an d the formation of vortices [331 can now be 
compared with experimental data |S1 • Needless to say that this dramatic progress on 
the experimental front has stimulated a wave of activity on both the theoretical and 
the numerical front. 

The properties of a BEC at temperatures T much smaller than the critical con- 
densation temperature T c j2HJ are usually well modeled by a nonlinear Schrodinger 
equation (NLSE), also called Gross-Pitaevskii equation (GPE) |29IK-tM| . for the macro- 
scopic wave function which incorporates the trap potential as well as the interactions 
among the atoms. The effect of the interactions is described by a mean field which 
leads to a nonlinear term in the GPE. The cases of repulsive and attractive inter- 
actions - which can both be realized in the experiment - correspond to defocusing 
and focusing nonlinearities in the GPE, respectively. The results obtained by solving 
the GPE showed excellent agreement with most of the experiments (for a review see 
Plll7j). In fact, up to now there have been very few experiments in ultra-cold dilute 
bosonic gases which could not be described properly by using theoretical methods 
based on the GPE j2H][2EJ- Thus developing efficient numerical methods for solving 
GPE is very important in numerical simulation of BEC. 



* Department of Computational Science, National University of Singapore, Singapore 117543 
(bao<Scz3.nus.edu.sg). URL: http://www.cz3.nus.edu.sg/~bao/ Research is supported by the Na- 
tional University of Singapore grant No. R-151-000-027-112. 

* Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA, 
(shen@math.purdue.edu). Research is supported by NSF DMS-0311915 and he acknowledges sup- 
port by the OAP (Fellow-Inbound) Programme supported by A*STAR and National University of 
Singapore. 



1 



2 



W. Bao and J. Shen 



Recently, a series of numerical studies are devoted to the numerical solution of 
time-independent GPE for finding the ground states and of time-dependent GPE for 
determining the dynamics of BECs. To compute ground states of BECs, Bao and 
Du |H] presented a continuous normalized gradient flow (CNGF) with diminishing 
energy, and discretized it by a backward Euler finite difference (BEFD) method; 
Bao and Tang P| proposed a method which can be used to compute the ground and 
excited states via directly minimizing the energy functional; Edwards and Burnett 20 
introduced a Runge-Kutta type method; other methods include an explicit imaginary- 
time algorithm in a directly inversion in the iterated subspace (DIIS) in [3l)| 
and a simple analytical type method in [ID]. To determine the dynamics of BECs, 
Bao et al. ^3 El presented a time-splitting spectral (TSSP) method, Ruprecht 
et al. used the Crank-Nicolson finite difference (CNFD) method, Cerimele et al. 
|14M15| proposed a particle-inspired scheme. 

In most experiments of BECs, the magnetic trap is with cylindrical symmetry. 
Thus, the 3-D GPE in Cartesian coordinate can be reduced to an effective 2-D prob- 
lem in cylindrical coordinate. In this case, both the TSSP |H3 EH IHj and CNFD 
[3"5] methods have serious drawbacks: (i) One needs to replace the original whole 
space by a truncated computational domain with an artificial (usually homogeneous 
Dirichlct boundary conditions are used) boundary condition. How to choose an ap- 
propriate bounded computational domain is a difficult task in practice: if it is too 
large, the computational resource is wasted; if it is too small, boundary effect will 
lead to wrong numerical solutions, (ii) The TSSP method is explicit and of spectral 
accuracy in space, but one needs to solve the original 3-D problem due to the pe- 
riodic/homogeneous Dirichlet boundary conditions required by Fourier/sine spectral 
method. Thus, the memory requirement is a big burden in this case. The CNFD 
discretizes the 2-D effective problem directly, but it is implicit and only second-order 
accurate in space. The aim of this paper is to develop a numerical method which 
enjoys advantages of both TSSP and CNFD. That is to say, the method is explicit 
and of spectral order accuracy in space, and discretizes the effective 2-D problem 
directly. We shall present such an efficient and accurate numerical method for dis- 
crctizing 3-D GPE with cylindrical symmetry by applying a time-splitting technique 
and constructing appropriately scaled Laguerre-Hermite basis functions. 

The paper is organized as follows. In Section [5J we present the Gross-Pitaevskii 
equation and its dimension reduction. In Section we present time-splitting Her- 
mite, Laguerre and Laguerre-Hermite spectral methods for 1-D GPE, 2-D GPE with 
radial symmetry and 3-D GPE with cylindrical symmetry, respectively. Extension 
of the time-splitting Laguerre-Hermite spectral method for vector Gross-Pitaevskii 
equations (VGPEs) for multi-component BEC is presented in Section 0] In Section 
numerical results for 1-D GPE, 2-D GPE with radial symmetry, 3-D GPE with 
cylindrical symmetry as well as 3-D VGPEs for multi-component BEC are reported 
to demonstrate the efficiency and accuracy of our new numerical methods. Some 
concluding remarks are given in Section [BJ 

2. The Gross-Pitaevskii equation (GPE). At temperatures T much smaller 
than the critical temperature T c [221, a BEC is well described by the macroscopic 
wave function tp = ip(pt,t) whose evolution is governed by a self-consistent, mean 
field nonlinear Schrodinger equation (NLSE) known as the Gross-Pitaevskii equation 

(gpe) using 

m tWgQ = _^v 2 V(x, t) + F(xMx, t) + NUoMx, t)\ 2 ^(x, t), (2.1) 
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where m is the atomic mass, h is the Planck constant, N is the number of atoms in the 
condensate, V^(x) is an external trapping potential. When a harmonic trap potential 
is considered, V(x) = ^ (w 2 x 2 + w 2 y 2 -f-w 2 z 2 ) with u x , lu v and uj z being the trap 
frequencies in x, y and z-direction, respectively. In most current BEC experiments, 
the traps are cylindrically symmetric, i.e. u> x = uj y . Uq — 4irh 2 a s /m describes 
the interaction between atoms in the condensate with the s-wave scattering length 
a s (positive for repulsive interaction and negative for attractive interaction). It is 
convenient to normalize the wave function by requiring 

|^(x,t)| 2 dx = 1. (2.2) 

2.1. Dimensionless GPE. In order to scale the Eq. 1(2. 1() under the normal- 
ization 1)2.2(1 . we introduce 

t = u} m t, x=— , -0(x,t) = a^/ 2 ip(x,t), with a = ^h/muj m , (2.3) 
ao 

where uj rn = min{w x , u v , u> z }, &o is the length of the harmonic oscillator ground 
state. In fact, we choose l/oj m and ao as the dimensionless time and length units, 
respectively. Plugging (|2.3|) into (|2.1|) . multiplying by l/moj^a^ 2 , and then removing 
all ~, we get the following dimensionless GPE under the normalization (|2.2|l in three 
dimension 

1 = *) + ^W^( x ' *) + I 3 W x > *)l^( x ' *)' ( 2 - 4 ) 

where 8 = = and 

V(x) = ~ (-f x x 2 + "fly 2 + j z z 2 ) , with j a = — (a = x, y, z). 

There are two extreme regimes of the interaction parameter 8: (1) 8 — o(l), the 
Eq. 1)2.4(1 describes a weakly interacting condensation; (2) 8 3> 1, it corresponds to a 
strongly interacting condensation or to the semiclassical regime. 

There are two typical extreme regimes between the trap frequencies: (1) 7^ = 1, 
7 y « 1 and j z 3> 1, it is a disk-shaped condensation; (2) 7^ 3> 1, 7^ ^> 1 and j z = 1, 
it is a cigar-shaped condensation. In these two cases, the 3-D GPE 1(2. 4() can be 
approximately reduced to a 2-D and 1-D equation respectively PHI El El as explained 
below. 

2.2. Reduction to lower dimension. Case I (disk-shaped condensation): 

uj x K,uj yi uj z ^uj x , 7 X = 1, 7 y 1, 7 Z > 1. 

Here, the 3-D GPE 12.4(1 can be reduced to 2-D GPE with x = (a;, y) by assuming that 
the time evolution does not cause excitations along the z-axis, since the excitations 
along the z-axis have large energy (of order Huj z ) compared to that along the x and 
y-axis with energies of order Hlu x . Thus, we may assume that the condensation wave 
function along the z-axis is always well described by the harmonic oscillator ground 
state wave function, and set 

iP = iP 2 (x,y,t)<p ho (z) with ho (z) = ( 7z A) 1/4 e-^ z2 ' 2 . (2.5) 
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Plugging (|2.5|) into l|2.4|) . multiplying by 4>^ (z) (where /* denotes the conjugate of a 
function /), integrating with respect to z over (—00,00), we get 

. a^t) = _ 1 ^ + 1 ^ + ^ 2 + c)ip2 + m ^ ( 2. 6) 



dt 



where 

(32= (3 



/OO I pOC / 

jL{z)dz = p S J^, C= J ^ ( 



l 2 z z 2 \^o(z)\ 2 



d(f) ho 



dz 



dz. 



Since this GPE is time-transverse invariant, we can replace ip2 - > e * 2 so that the 
constant C in the trap potential disappears, and we obtain the 2-D effective GPE: 



dt 



(2.7) 



Note that the observables, e.g. the position density |V>| 2 , are not affected by dropping 

the constant C in (|2.tj|) . 

Case II (cigar-shaped condensation): 



Wl » My > LJ 



7^ > 1, 7y > 1> lz = 1- 



Here, the 3-D GPE (12.411 can be reduced to a f-D GPE with x = z. Similarly as in 
the 2-D case, we can derive the following f-D GPE |5ffl HTH ^ : 

1 ^aT~ = ~l^ z > *) + ^yV*. *) + AM*. 0, (2-8) 

where ft = (3^%7^/2tt. 

The 3-D GPE 2-D GPE J23) and f-D GPE can be written in a unified 

form: 



■0(x,O) = V>o(x), 



x e 



(2.9) 
(2-10) 



with 



y / 7^/2 7T, f 7^ 2 A 

fti=/^ VW^- ^ d (x) = ^ ( 7 2 x 2 + 7 2 y 2 ) A 

1 f, I (7 2 ^ 2 + 7|y 2 + llz 2 ) /2, 



d = l, 
d = 2, 
d = 3, 



where 7 X > 0, 7a > and 72 > are constants. The normalization condition for (|2.9|l 



W(V0 = lhK-,*)l| 2 = / |^(x,t)| 2 dx= / |^o(x)| 2 dx=f 



(2-11) 



3. Fourth-order time-splitting Laguerre-Hermite pseudospectral method.| 

In this section we present a fourth-order time-splitting Laguerre-Hermite pseudospec- 
tral method for the problem l|2.9|) - (|2.f U|l in 3-D with cylindrical symmetry. As 
preparatory steps we begin by introducing the fourth-order time-splitting method 
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and applying it with Hermite pseudospectral method for 1-D GPE and with Laguerre 
pseudospectral method for 2-D GPE with radial symmetry, respectively. 
Consider a general evolution equation 

iu t = f(u) = Au + Bu (3.1) 

where f(u) is a nonlinear operator and the splitting f(u) = Au + Bu can be quite 
abitrary, in particular, A and B do not need to commute. For a given time step 
At > 0, let t n = n At, n — 0,1,2,... and u n be the approximation of u{t n ). A 
fourth-order symplectic time integrator (cf. for (|3.1[) is as follows: 

u (l) = e ^2 Wl AAt u n. 
u (2) = e -t2 W2 BAt u (l) 



„(3) = e -i2w 3 AAt u {2) 
u (4) = e -i2 Wi BAt u (3) 
u (5) = e -i2 W3 AAt U (A) 
u (6) = e -i2 W2 BAt u (5) 
u n+l = e -i2 wl AAt u (6) 



(3.2) 



where 



wi = 0.33780 17979 89914 40851, w 2 = 0.67560 35959 79828 81702, 
w 3 = -0.08780 17979 89914 40851, w A = -0.85120 71979 59657 63405. 

We now rewrite the GPE l|2.9|l in the form of (|3.1|l with 



(3.3) 



Aip=f3 d ^(x,t)| 2 ^(x,i), B^ = -iv 2 ^(x,t)+V r tI (x)V(x,t). (3.4) 

Thus, the key for an efficient implementation of (|3.2|l is to solve efficiently the following 
two subproblems: 



and 



* = ^(x, t) = -± VV(x, t) + ^(x)^(x), x e M d , 

lim V( x > *) = 0. 



(3.6) 



The decaying condition in l|3.6|) is necessary for satisfying the normalization 12.11fl . 

Multiplying (|3.5() by 4>{x., t), we find that the ordinary differential equation (|3.5|l 
leaves |^(x,i)| invariant in t. Hence, for t>t s (t s is any given time), (|3.5|l becomes 

1 &t =(3d l^*^)! 2 ^*^ xeMd ( 3 - 7 ) 

which can be integrated exactly, i.e., 

ii(x,t)=e- ll3d ^' t ^ 2( > t - t ^iP(x,t s ), t>t s , xeR d . (3.8) 
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Thus, it remains to find an efficient and accurate scheme for (|3.6[) . We shall con- 
struct below suitable spectral basis functions which are eigenfunctions of B so that 
e -iBAt^ can k e exactly evaluated (which is necessary for the final scheme to be time 
reversible and time transverse invariant) . Hence, the only time discretization error of 
the corresponding time splitting method (|3.2|l is the splitting error, which is fourth 
order in At. Furthermore, the scheme is explicit, time reversible and time transverse 
invariant, and as we shall show below, it is also unconditionally stable. 

3.1. Hermite pseudospectral method for the 1-D GPE. In the 1-D case, 
Eq. (|3.fi[) collapses to 

dvb(zA) , . ld 2 ib(z,t) ^ 2 z 2 , , , m . . 

at = t] = ~ 2 dz 2 + 2 t); ze ' ( } 

lim ip[z,t) =0,t>0, (3.10) 

\z\ — > + oo 

with the normalization (|2.11|) 

U(;t)\\ 2 = / \^(z,t)\ 2 dz = Wo(z)\ 2 dz = l. (3.11) 



Since the above equation is posed on the whole line, it is natural to consider Hermite 
functions which have been successfully applied to other equations (cf. [221 EH)- Al- 
though the standard Hermite functions could be used as basis functions here, they 
are not the most appropriate. Below, we construct properly scaled Hermite functions 
which are eigenfunctions of B. 

Let Hi(z) (I = 0,1, ... , N) be the standard Hermite polynomials satisfying 

H'i'(z)-2zHl(z)+2lHi(z) = 0, zel, / > 0, (3.12) 

/oo 

Hi(z)H n (z)e- z2 dz = V?F 2 z Z! S tn , l,n>0, (3.13) 
-OO 

where Si n is the Kronecker delta. We define the scaled Hermite function 

h l (z)^e^ z2 / 2 H l (^z)/V¥v.(7r/ lz ) 1 /\ zel. (3.14) 

Plugging [j3~TlT into [ETTlj) and ll3~T3|) . we find that 

1 ^ 2 z 2 21 + 1 

--hnz) + 2 Yh l (z)=tfh l (z), zeR, M f = _Z^ 7 „ Z>0, (3.15) 

/°° 1 2 

-.H l {z)H n {z)e' z dz = 5 ln , l,n>0. (3.16) 
-oo Vir2 l l\2 n nl 

Hence, {hi} are eigenfunctions of B defined in H3.9J) . 

For a fixed N, let Xn = spanj/i; : / = 0, 1, • • • , N}. The Hermite-spectral method 
for l|3.9|l is to find ipN(z,t) S Xn, i.e., 

N 

ip N (z,t) = Y,Mt) h{z), zel. (3.17) 

1=0 

such that 

■ dip N {z,t) l d 2 4> N {z,t) -ilz 2 

i — = B<ip N (z,t) = — 2 + —^-ifj N (z,t), z£l. (3.18) 
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Note that lim| 2 |^ +00 h t (z) — (cf. 0U1) so the decaying condition lini| z |^ +00 iPn{z> t) 
= is automatically satisfied. 

Plugging (|3T7|l into l|3TT8)l . thanks to I^TH)} and foT^)l . we find 



, #P = Mf ^ (t) = ^[ + l 7z ^ i = 0, 1, • • ■ , JV. (3.19) 



Hence, the solution for (|3.18l) is given by 



JY 



lM*>*) = e- iB ^-^ip N (z,t s ) =^ e - i "?(*-'.) ^(tMz), t>t s . (3.20) 



J=0 



Let {^fcjfcLo be tne Hermite-Gauss points (cf. @U1E21)i i- e - {zk}k=o are * ne JV+1 
roots of the polynomial H N+ i(z). Let ^>J? be the approximation of tp{z k ,t n ) and i\> n 
be the solution vector with components ip k l . Then, the fourth-order time-splitting 
Hermite-pseudospectral (TSHP4) method for 1-D GPE (|2.9|l is given by 



4 



4 



4 



-i2wi At /3i|Vfc| 2 „/,n 



N 



= J2 e ~ i2W2 Mf At W (1) ) ; M*k), 

1=0 

_ e -i2w 3 AtftlV-i 2 '! 2 ^,(2) ) 
JV 

= J2e~ i2w ^' At tyW), h{z k ), fc = 0,l 

_ -i2iu 3 At /9i | v£" 



1=0 



? 5 * * *? ? 



JV, 



(3.21) 



where u>j, i = 1, 2, 3,4 are given in (|3.3|) . and {I 7 /}, the coefficients of scaled Hermite 
expansion of U(z), can be computed from the discrete scaled Hermite transform: 



JV 



Ui = Y, ^kU(z k )hi(z k ), 1 = 0,1,. ..,N. 



(3.22) 



fc=0 



In the above, z k and ui k are the scaled Hcrmitc-Gauss points and weights, respectively, 
which are defined by 



Zk 



"k 



< fc < JV, 



(3.23) 



where {w^}^ are the weights associated with the Hermite-Gauss quadrature (cf. 
|22p satisfying 



JV 



Hi(zk) H„(z k ) 



k=0 



Sin, l,n = 0, 1 



,...,JV, 



(3.24) 
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and we derive from 1|3.14|) that 

JV JV 

^ luI hi(z k ) h m (z k ) = ^ ^k e7 ' k l^/Tz hi {z k f^u) h m {zk/y/%, 

Hi(z k ) H n (z k ) 



k=0 



k=0 
JV 

E 

k=0 



= Sin, 0<l,n<N. (3.25) 



Note that the computation of {uJ k } from l|3.23|) is not a stable process for very large 
N. However, one can compute {^ k } in a stable way as suggested in the Appendix of 

m 

Thus, the memory requirement of this scheme is O(N) and the computational 
cost per time step is a small multiple of N 2 . As for the stability of the TSHP4, we 
have the following 

Lemma 3.1. The time- splitting Hermite-pseudospectral (TSHP4-) method \S. 21)) 
is unconditionally stable. More precisely, we have 



JV 



M 



wwi = E«i 2 = E*°(^)i 2 = iiV'oi 



n = 0,1, 



(3.26) 



k=0 



k=0 



Proof. From noting $£22$ and we obtain 



JV 



JV 



= E«i 2 = E 



} z L-»2u>i At ft | i,[ 6) | 2 ^(6) 



fe=0 
JV 



k=0 
N 



E^i€ } i 2 =E^ 



fc=0 
JV JV 

EE- 

Z=0 m=0 
JV JV 

EE- 

1=0 m=0 
JV 



fc=0 



-z2u!2 At 



-z2u!2 l±i At 



JV 



JV 

E< 



-ilw-2 \i\ Ai 



(^ (5) )j M**) 



-&2u>2 jU™ At 



(W> (5) )J* 



JV 



^uj k hi(z k )h m (z k ) 



,fc=0 



(^(5)) e ^ 2 „5. At ((^(5)) )* ^ 



J=0 
JV JV 



Z=0 



JV 



Ei(^ (5) )ji 2 = E E ^^ (5, wm«) 



Jc=0 



JV 



fe=0 m=0 
JV JV 

fc=0 m=0 
JV 

E«2l^ B) (ft)l a = ll^ (5, llp 



^LU^hi(z k )hi(z m ) 



.1=0 



km 



k=0 

Similarly, we have 



(3.27) 



(3.28) 
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Thus the equality l|3.26[) can be obtained from (|3.28|) by induction. □ 

Remark 3.1. Extension of TSHP4 method to 2-D GPE without radial 

symmetry and 3-D GPE without cylindrical symmetry is straightforward by using 
tensor product of scaled Hermite functions. 

3.2. Laguerre pseudospectral method for 2-D GPE with radial sym- 
metry. In the 2-D case with radial symmetry, i.e. d = 2 and ^ x = j y in (|2.9|l . and 
tpo(x,y) = tpo(r) m 1|2.10[1 with r = \J 1 x 2 + y 2 , we can write the solution of i|2.9|) . 
(|2.10|) as 4>(x,y,t) = ip(r,t). Therefore, Eq. (|3.6j) collapses to 

dMr,t) . . 1 d ( dMr,t)\ j 2 r 2 . . , 

1 dT~ = mr ' t] = -Yrd~r { r ^r^J + ° < T < °°' (3 ' 29) 

lim Mr, f) = 0, t>0, (3.30) 

r — >oo 

where 7 r = j x = j y . The normalization (|2.11|1 collapses to 

||^(-,*)|| 2 = 2tt / \ip(r,t)\ 2 r dr = 2ir \i/j (r)\ 2 r dr = 1. (3.31) 



Note that it can be shown, similarly as for the Poisson equation in a 2-D disk (cf. 
GUI), that the problem (jS^-ESOJ admits a unique solution without any condition 
at the pole r = 0. 

Since (|3.29() is posed on a semi-infinite interval, it is natural to consider Laguerre 
functions which have been successfully used for other problems in semi-infinite inter- 
vals (cf. (221 OH] ) - Again, the standard Laguerre functions, although usable, are not 
the most appropriate for this problem. Below, we construct properly scaled Laguerre 
functions which are eigenfunctions of B. 

Let L m (r) (m = 0, 1, . . . , M) be the Laguerre polynomials of degree m satisfying 

rL'^(r) + {l-r)L' m (r)+mL m (r) = 0, m = 0,l,..., (3.32) 
e~ r L m (r) L n (r) dr = 5 mn , m,n = 0, 1, . . . . (3.33) 

We define the scaled Laguerre functions L m by 

L m (r) = ^ e-^ r2/2 L m ( lr r 2 ), < r < oo. (3.34) 

Note that lim| r |^ +oc L m (r) = (cf. 140] ) hence, limi r u +KI ipM(f, t) — is automati- 
cally satisfied. 

Plugging H3.34J1 into 13.32J1 and 1)3.33(1 . a simple computation shows 
-~ + \iyLm{r) = H r m L m (r), tf m = lr {2m + 1), m > 0, (3.35) 

y* oo />oc 

2tt L m (r)L n (r)r dr = / e~ r L m (r)L n (r) dr — 5 mn , m,n>0. (3.36) 
Jo Jo 

Hence, {L m } are eigenfunctions of B defined in 13.29fl . 

For a fixed M, let Ym — span{L m : m = 0, 1, • • • , M}. The Laguerre-spectral 
method for l|3.9|l is to find ipM{r,t) € Ym, i.e., 

M 

^m{t, t)=J2 L ™W' < r < oo (3.37) 

m=0 
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such that 
. dtp M (r,t) 



Or 



1 d ( 9VmM)\ , l 2 r r 2 



+ J: Y-ii>M{r,t), 0<r<oo. 



Plugging ljX?7jl into [j3~351l . thanks to (|3~35|) and 133(3, wc find 
. dip m (t) 



dt 



\i r m i> m {t) = 7 z (2m + l)V'm(i), m = 0, 1, 



(3.38) 
(3.39) 



Hence, the solution for (|3.38|l is given by 



M 



ipM{r,t) = e 



- e -iB(t-t, 



] ipM(r,t s ) = ^2 



^m(t s )L m (r), t>t s . (3.40) 



m=0 



Let {rj}fL be the Laguerre-Gauss-Radau points (cf. jUJ), i.e. they are the M+l 
roots of the polynomial rL' M+1 (r). Let ipj be the approximation of ip(rj,t n ) and ip n 
be the solution vector with components ip™. Then, the fourth-order time-splitting 
Laguerre-pseudospectral (TSLP4) method for 2-D GPE l|2.9|l with radial symmetry 
is given by 



(i 



(2 



(3 



(4 



(5 



(6 



,1/ 



i=0 



-i2ro 3 At/3 2 |^ 2) | 2 ^(2) 



= E e ^ 2W4 * At (^ (3) )j L '( r i)' 3 = 0, 1, • ■ • , M, 



1=0 



e -i2w 3 At/3 2 |</>< 4) | 2 ^(4)^ 



2=0 



(3.41) 



^n+l = e -i2iu 1 At /9 2 |Vf 'I 2 ^( 6 ). 

where t/i , the coefficients of scaled Laguerre expansion of U(r) can be computed from 
the discrete scaled Laguerre transform: 

M 

Ui = ^2^U(r j )L l (r j ), l = Q,l,...,M. (3.42) 

In the above, Vj and w| are the scaled Laguerre-Gauss-Radau points and weights, 
respectively, which are defined by 



r _ Q r r 3 _ 



i = 0,l,...,M, 



(3.43) 



where {uJj}fL are the weights associated to the Laguerre-Gauss quadrature [22] sat- 



isfying 



M 



}^w'-L m (fj)L n (fj) = 5 nm , n,m = 0,1, ...,M, 
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and we derive from l|3.34|) that 



M 



M 



3=0 



3=0 
M 



= UjLm(rj)L n {rj) = S nm , n,m = 0,l,...,M. (3.44) 

3=0 

As in the Hermite case, the computation of {^J} from l|3.43[l is not a stable process 
for very large N. However, one can compute {wj} in a stable way as suggested in the 
Appendix of |3"%] . 

The memory requirement of this scheme is O(M) and the computational cost per 
time step is a small multiple of M 2 . As for the stability of the TSLP4, we have the 
following 

Lemma 3.2. The time- splitting Laguerre-pseudospectral (TSLP4-) method \3.41\j 
is unconditionally stable. More precisely, we have 

M M 

wrwi = $>Mi 2 = £^rM^)| 2 = n > 0. 

3=0 3=0 



Proof. Using (|3.44|) , the proof is essentially the same as in Lemma 13.11 for time- 
splitting Hermite-pseudospectral method. □ 

3.3. Laguerre-Hermite pseudospectral method for 3-D GPE with cylin- 
drical symmetry. In the 3-D case with cylindrical symmetry, i.e. d — 3 and j x = j y 
in H2.9J) . and ipo(x,y,z ) = ipo(r,z) in (CTty . the solution of $n fy -$n fy with d = 3 
satisfies i/j(x,y, z,i) = ip(r,z,t). Therefore, Eq. 1)3. 6|l becomes to 



» m = B^(r,z,t) = -- 



1 d ( dip\ d 2 ip 
r dr \ dr J dz 2 



< r < oo, — oo < z < oo, 
lim ip(r, z, t) = 0, lim ^>(r, z,t) = 0, t > 0, 



(3.45) 
(3.46) 



where 7 r = j x — j v . The normalization (|2.11[1 becomes 
U(-,t)\\ 2 = 2n 



oo poo 



oo poo 



\tp(r, z,t)\ 2 r drdz = 2tt / / \ipo(r, z)\ 2 r drdz = 1. 

'0 J— oo JO J -oo 

(3.47) 

We are now in position to present our Laguerre-Hermite pseudospectral method for 

{EE} . 

Using the same notations as in previous subsections, we derive from i|3.15[) and 
that 



19/9 



d 2 



r dr V dr I dz 



(L m (r) h{z)) + ~ ( 7r V + ~t 2 z z 2 ) (L m (r) &,(*)) 



1 d fJLm{r ^ + ^y Lm{ r) 



2r dr 



dr 



hi(z) + 



l d 2 hi{z) 1 22 ' 
'2^^ + 2 7 ^ /l ' (z) 



/J^Lmi^hiiz) + fifh l (z)L m (r) = [fj, r m + nf)L m (r)hi(z). 



L m (r) 
(3.48) 
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Hence, {L m (r)hi(z)} are eigenfunctions of B denned in 

For a fixed pair (M,N), let Xmn = span{L m (r)/ii(z) : m — 0, 1, • • • , M, I = 
0, 1, • ■ • , N}. The Laguerre-Hermite spectral method for (|3.45[) is to find iPmn(t, z, t) 
S Xmn, i-e., 



M N 



ipMN(r,z,t) = ^ /^ml(t) L rn (r) hi(z) 



(3.49) 



m=0 1=0 



such that 

. dtpMNjr, z,t) 

1 at 



~ld_ 

r dr 



1 d ( dipMN 
dr 



+ 



dz 2 



+ 2 hr 7 * 2 +ll z2 ) ^MN- 



Plugging iPH9"|) into [|3"3Ull . thanks to (|3~4^|) . we find that 

. dijj m i(t) 



dt 



{ l f m + Ht)^ ml {t), m = 0,l,--- ,M, 1 = 0,1,. 



,iV. 



(3.50) 



(3.51) 



Hence, the solution for (|3.5UI) is given by 

ipMN{r,z,t) = e' lB(t - ts) iljMN{r,z,t s ) 



M N 



+Mf)(*-t s ) 



4> m i(t s )L m (r) h(z), t>t s . (3.52) 



m=0 1=0 



Let ip™ k be the approximation of ip{rj, z kl t n ) and ip n be the solution vector with 
components The fourth-order time-splitting Laguerre-Hermite-pseudospectral 

(TSLHP4) method for 3-D GPE (|2.9|) with cylindrical symmetry is given by 



_ -i2u;i At ,931V?*! 



A/ N 



m=0 1=0 



e -i2u, 3 At/3 3 |</>< 2) | 2 ^(2) 



EE' 



e -i2. 4 A t (^ + ^) ( ^ (3))mj Lm ( rj)/li ( Zfe)) 

e -^ 3 At ftl^fl 2 ^ j = 0, 1, . . . , M, fc = 0, 1 . . . , AT, 

M N 

J- ^y^At^+tf) ( ^5) )mj jLm(rj)/li(zfc)j 
m=0 2=0 

e~ 



= p- l2t01 Atfth^T ,/.(6). 



(3.53) 



-jfc ' rjk ' 

where C/ m /, the coefficients of scaled Laguerre-Hermite expansion of U(r, z) are com- 
puted by the discrete scaled Laguerre-Hermite transform 

M N 

Uml = /]/,Vj u k U(rj,z k ) L m (rj)hi(z k ), m = 0,l,...,M, k = 0, 1, . . . , N. 



j=0 k=0 



(3.54) 
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The memory requirement of this scheme is O(MN) and the computational cost per 
time step is 0(max(M 2 N, N 2 M)). As for the stability of the TSLHP4, we have the 
following 

Lemma 3.3. The time- splitting Laguerre-Hermite pseudospectral (TSLHP4 ) method^ 
13. 5S\) is unconditionally stable. More precisely, we have 

M N M N 

j=0 k=0 3=0 k=0 



Proof. Using l|3.25l) and 13.44|l . the proof is essentially the same as in Lemma l3~Tl 
for time-splitting Hermite-pseudospectral method. □ 

4. Extension to multi-component BECs. The time-splitting Laguerre-Hermite| 

pseudospectral method, introduced above for the 3-D GPE with cylindrical symmetry, 
can be extended to vector Gross-Pitaevskii equations (VGPEs) for multi-component 
BECs 5 . For simplicity, we only present the detailed method for the dynamics of 
two-component BECs. Consider the dimensionless VGPEs with an external driven 
filed (cf. 0) 



dip{r,z,t) 
dt 



ld_ 

r dr 



dr 



d 2 ip 
dz 2 



+ (/3nH 2 + /3 12 |0| 2 ) V + y/N$/N°f(t)<l>, 



d(j){r,z,t) 
dt 



ld_ 

r dr 



dr 



dz 



+ ~ 2 h> 2 



,0\2\ 



(4.1) 



(4.2) 



+ (/32i|Vf + P2M 2 ) + \jN°/N°f(t)i/>, < r < 00, z € R, 
lim ip(r,z,t) = lim cj>(r, z, t) = 0, lim tf)(r,z,t) = lim (f>(r, z,t)=0, (4.3) 

r^oo r^oo \z\^oo |z|^oo 

i/j(r,z,0) = ipa(r,z), 



(r, z,0) = <t> {r, z); 



(4.4) 



where z® and N® (j = 1,2) are the center of trapping potential along z-axis and 
the number of atoms of the jth component, respectively; j r — u) r /u> m , 7 Z = uj z /uj m 
with io r , lo z and oj m are the radial, axial and reference frequencies, respectively; (3ji — 
AirajiN® /ao — 1,2) with the s-wave scattering length dji = aij between the jth 
and Zth component and a<j = ^Jh/mL0 m ] and fit) = cos{uidt / u>m) / u m with f2 and 
ojd being the amplitude and frequency of the external driven field. The wave functions 
are normalized as 



OO fOO 



= 2?r / / \tp {r,z)\ 2 r drdz = 1, 

JO J — oo 

j* OO />00 

||0o|| 2 = 2W / \Mr,z)\ 2 rdrdz = l. 
Jo J -00 

It is easy to show (cf. 5 ) that the total number of atoms is conserved 



(4.5) 



N° 1 M;t)\\ 2 + N°H(;t)\\ 2 = 2n 



00 poo 



(\ip(r,z,t)\ 2 -I- \(f)(r,z,t)\ 2 ) r drdz 



N^UoW 2 + N°\\M 2 = 



N° 2 . 



(4.6) 
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Unlike the time-splitting Laguerre-Hermite method for 3-D GPE (|2.9f) with cylindrical 
symmetry, here we have to split the VGPEs l|4. 111 - 14. 2f) into three sub-systems. For 
example, for a fist-order splitting scheme, we first solve 



. dip{r,z,t) 
% dt 
■ dc/>(r,z,t) 



id_ 

r dr 
ld_ 

r dr 



(hp 
dr 
dfi 
dr 



d 2 *p 
dz 2 
d 2 f 
dz 2 



2 US 



2 2, 2 2\ i 



(4.7) 
(4.8) 



for the time step of length At, followed by solving 



t d^(r, z,t) = l 7 2 z?(z o _ 2z) ^ + + ^ M2) ^ 



dt 

94>{r, z, t) 
dt 



1 



7 2 z z°(z° 2 - 2z)<j> + (/? 21 |Vf + P2M 2 ) & 



(4.9) 
(4.10) 



for the same time step, and then by solving 



(4.11) 
(4.12) 



for the same time step. 



The nonlinear ODE system (|4.9|) - (|4.1U|) leaves \ip(r,z,t)\ and \<fi(r,z,t)\ invariant 
in t and thus can be integrated exactly The linear ODE system 14.llfl - H4.12fl 
can also be integrated exactly by applying a matrix diagonalization technique (cf. 
[1]). As is shown above, (|4.7fl - (|4.8f) can be discretized in space by Laguerre-Hermite 
pseudospectral method and integrated in time exactly. 



Let ipj'k an d 4*jk b e the approximations of ij)(rj, Zk,t n ) an d 4>{rj,Zk,t n ), respec- 
tively, and ip n and 4> n be the solution vectors with components tp™ k and respec- 
tively. Although it is not clear how to construct a fourth-order time splitting schemes 
with three sub-systems, a second-order scheme can be easily constructed using the 
Strang splitting (cf. |39|). More precisely, from time t = t n to t = t n +i, we proceed 



Laguerre-Hermite pseudospectral method for BEC 



15 



as follows: 

M N 

',(1 



m=0 i=0 
M AT 



m=0 Z=0 

,/.(2) = p -i[7^?(2?-2^)/2+(/3ii|V'jt ) | 2 +^ 2 |4 1 fe ) | 2 )]At/2 / (l) 
,(2) _ - l [ 7 f z «(4_2 Zfc )/2+(fe 1 |V^ ) | 2 +/3 22 |0< 1 fc ) | 2 )]At/2 ,(1) 

1$> = cos( ff (i n+1 ,t„))^ 2 fe ) -ism(s(t n+1 ,t n ))yjN$/N?<l$, 
<f>$ = -i sm(g(t n+1 ,t n ))yj N°/N°il>™ + cos( 5 (i n+1 , t„ ))0$ , 

,/,( 4 ) - p - l [7f2 1 (^ 1 -2)/2+(/3 11 |vf fc ) | 2 +/3 12 |0f fc ) | 2 )]At/2 ,(3) 

(4) = e -<[^.S(.S-a)/^(fri|^>| a +/fe a |*W| a )]At/2 (8) i o < j < M, < fc < TV, 

M JV 

= E ^e-^+^)At/2 (^) m; Lm{rj)hl{zk)i 

rn=0 1=0 
M N 

£ ^ e -^ m+ Mf)At/2 Lm(r . )/li( ^ )) (4.13) 



m=0 i=0 



where 



g(t,t n ) = f(s) ds = Qui d [shx{udt/uj m ) - sm(cj d t„/uj m )} . 

Note that the only time discretization error of this scheme is the splitting error, which 
is of second order in At. The scheme is explicit, spectral accurate in space and second 
order accurate in time. The memory requirement of this method is O(MN) and the 
computational cost per time step is 0(max(M 2 N, MN 2 )). As for the stability, we 
can prove as in [5] the following lemma, which shows that the total number of atoms 
is conserved in the discretized level. 

Lemma 4.1. The time- splitting Laguerre-Hermite-pseudospectral method \4-liH\j 
for multi- component BEC is unconditionally stable. More precisely, we have 



5. Numerical results. We now present some numerical results by using the 
numerical methods introduced in section [3] To quantify the numerical results, we 
define the condensate width along r- and z-axis as 

a l= a 2 \ijj(x,t)\ dx, a = x,y,z, a 2 = a 2 c + a 1 . 

Example 1. The 1-D Gross-Pitaevskii equation: We choose d = 1, j z = 2, 
Pi = 50 in i|2.9|) . The initial data ^o(z) is chosen as the ground state of the 1-D GPE 
(|2.9|) with d = 1, 7 Z = 1 and f3% — 50 [HUH!- This corresponds to an experimental 
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a). t b). t 

Figure 1: Evolution of central density and condensate width in Example 1. ' — ': 'exact 
solutions' obtained by TSSP ^U] with 1025 grid points over an interval [—12, 12]; '+ 
+ + ': Numerical results by TSHP4 (|3.21|) with 31 grid points on the whole z-axis. 
a). Central density \ip(0, t)\ 2 ; b). condensate width a z . 

setup where initially the condensate is assumed to be in its ground state, and the 
trap frequency is double at t = 0. We solve this problem by using (|3.21|) with N = 31 
and time step k = 0.001. Figure 1 plots the condensate width and central density 
\tp(0, t)\ 2 as functions of time. Our numerical experiments also show that the scheme 
(I3.21|) with N — 31 gives similar numerical results as the TSSP method |101 ITT] for 
this example with 129 grid points over the interval [—12, 12] and time step k = 0.001. 

Example 2. The 2-D Gross-Pitaevskii equation with radial symmetry: we choose 
d = 2, 7 r = j x — 7 y = 2, 02 = 50 in (|2.9(l . The initial data ipo(r) is chosen as the 
ground state of the 2-D GPE with d = 2, j r = j x = j y = 1 and = 50 

[HI IS]' Again this corresponds to an experimental setup where initially the condensate 
is assumed to be in its ground state, and the trap frequency is doubled at t = 0. We 
solve this problem by using l|3.41|l with M — 30 and time step k = 0.001. Figure 2 
plots the condensate width and central density |-i/>(0,i)| 2 as functions of time. Our 
numerical experiments also show that the scheme (|3.41|l with M — 30 gives similar 
numerical results as the TSSP method [lOHllj for this example with 129 2 grid points 
over the box [—8, 8] 2 and time step k = 0.001. 

Example 3. The 3-D Gross-Pitaevskii equation with cylindrical symmetry: we 
choose d = 3, j r — -f x = j y = 4, j z = 1 and (}% — 100 in l|2.9|) . The initial data ipo( r , z ) 
is chosen as the ground state of the 3-D GPE (|2.9() with d = 3, 7 r = j x = ~f y = 1, 
72 = 4 and /?3 = 100 jHIE]- This corresponds to an experimental setup where initially 
the condensate is assumed to be in its ground state, and at t = 0, we increase the 
radial frequency four times and decrease the axial frequency to its quarter. We solve 
this problem by using l|3.53|l with M — 60 and N — 61 and time step k = 0.001. 
Figure 3 plots the condensate widths and central density 1^(0, 0, t)\ 2 as functions of 
time. 

The numerical results for these three examples clearly indicated that our new 
methods are very efficient and accurate. 

Example 4. The 3-D vector Gross-Pitaevskii equations with cylindrical symme- 
try for two-component BECs: we take, in (|4.1|l and l|4.2|) . m — 1.44 x 10 -25 [kg], a\ 2 = 
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a). t b). t 

Figure 2: Evolution of central density and condensate width in Example 2. ' — ': 
'exact solutions' obtained by TSSP JU| with 1029 2 grid points over a box [— 8,8] 2 ; 
'+ + + ': Numerical results by TSLP4 i|3.41[) with 30 grid points on the semi- infinite 
interval [0,oo). a). Central density \ip(0, i)| 2 ; b). condensate width a r . 




O 2 4 6 S 

t 

Figure 3: Evolution of central density and condensate width in Example 3 by the 
TSLHP4 (ET33II . 

0-21 = 55. 3A = 5.53 [nm], an = 1.03ai2 = 5.6959 [nm], a 2 2 = 0.97ai 2 = 5.3641 [nm], 
lu z = 47 x 2tt [1/s], uj m = uj r = lu z /V8, N% = 7V 2 ° = 500,000, Q = 65 x 2tt [1/s], 
u>d = 6.5 x 2tt [1/s]. A simple computation shows ao = 0.2643 x 10 -5 [m], (3n = 
0.02708165^°, (3 12 = 0.02629286^°, P21 = 0.02629286^°, (3 2 2 = 0.02 5 5 4 077V 2 °. The 
initial data ipo(r, z) and (f>o(r,z) are chosen as the ground state of the 3-D VGPEs 
(|4.1|) and (|4.2() . and we set f(t) = 5 . We solve this problem by using (|4.13|l with 
M = 100 and N — 201 and time step k = 0.00025. Figure 4 displays the time evolu- 
tion of the density functions for the two components with different trapping centers. 
The results are similar as those obtained in |Sj by a TSSP method with a much refined 
grid. 

^From Fig. 4, we can see that the general form of time evolution on the number 
of particles in the two components is similar for different distances between the two 
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Figure 4: Time evolution of the density functions for the two-component BECs in 
Example 4. a). z\ = z% = 0, b). z\ = -z% = 0.15, c). z\ = -z% = 0.4. 
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trapping potential centers. When z\ = — 0, the number of particles in the second 
component, i.e. iV^H^H 2 , decreases, reaching its bottom, oscillates and then attains 
its maximum at around t = 5.2. The number of particles in the second component 
at its maximum is approximately 55% bigger than its initial value at time t — 0. (cf. 
Fig. 4a). The pattern for 7V°||V>|| 2 is exactly the opposite of that for N° \\(f>\\ 2 (cf. Fig. 
4a). This antisymmetry is due to the fact that the total number of particles in the 
two components are conserved. When z® — z% > becomes larger, i.e. initially the 
density functions for the two components are separated further, the earlier the number 
of particles attains its absolute peak (cf. Fig. 4b&c), the smaller the maximum of 
the peak is. In fact, when z® — z% = 0.3 (resp. 0.8), at around t = 3.4 (resp. 
2.05), the number of particles in the first component attains its maximum which is 
approximately 52% (resp. 38.5%) bigger than its initial value at time t = 0. 

6. Concluding remarks. We developed a new efficient fourth-order time-splitting| 

Laguerre-Hermite pseudospectral method for 3-D Gross-Pitaevskii equation with cylin-| 
drical symmetry for Bose-Einstein condensates. The new method takes advantage of 
the cylindrical symmetry so only an effective 2-D problem is solved numerically. The 
new method is based on appropriately scaled Laguerre-Hermite functions and a fourth- 
order symplectic integrator. Hence, it is spectrally accurate in space, fourth-order 
accurate in time, explicit, unconditionally stable, time reversible and time transverse 
invariant. 

When Compared with the time-splitting sine-spectral method in |101 1111 [5] and 
the Crank-Nicolson finite difference method in |35l 12*] . the new method enjoys two 
important advantages: (i) there is no need to truncate the original whole space into 
a bounded computational domain for which an artificial boundary condition (which 
often erodes the accuracy) is needed; (ii) it solves an effective 2-D problem instead 
of the original 3-D equations. Thus, the new method is very accurate and efficient, 
particularly in term of memory requirement. Therefore, it is extremely suitable for 
3-D GPE with cylindrical symmetry which is the most frequently setup in BEC exper- 
iments. We plan to apply this powerful numerical method to study physically more 
complex systems like multi-component BECs, vortex states and dynamics in BECs. 
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